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Modeling and interpolation of the 
ambient magnetic fieid by Gaussian 
processes 
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Abstract 

Anomalies in the ambient magnetic field can be used as features in indoor positioning and navigation. By using 
Maxwell’s equations, we derive and present a Bayesian non-parametric probabilistic modeling approach for 
interpolation and extrapolation of the magnetic field. We model the magnetic field components jointly by imposing 
a Gaussian process (GP) prior on the latent scalar potential of the magnetic field. By rewriting the GP model in 
terms of a Hilbert space representation, we circumvent the computational pitfalls associated with GP modeling 
and provide a computationally efficient and physically justified modeling tool for the ambient magnetic field. The 
model allows for sequential updating of the estimate and time-dependent changes in the magnetic field. The 
model is shown to work well in practice in different applications: we demonstrate mapping of the magnetic field 
both with an inexpensive Raspberry Pi powered robot and on foot using a standard smartphone. 
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1 Introduction 

Ever since the invention of the compass, the direction 
of the Earth magnetic field has been used for navigation 
purposes by mankind. Moreover, animal behaviour has 
been influenced both by the direction and by the 
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intensity of the Earth magnetic field (see, e.g., Kirschvink 
and Gould 1981). Magnetic materials, however, cause 
anomalies in the ambient magnetic field. These materials 
are present on a large scale in indoor environments, for 
instance in the structures of buildings and large furniture. 
In this paper, we focus on building maps of the magnetic 
field. These maps are constructed by interpolating three- 
dimensional magnetic field measurements obtained using 
magnetometers. These sensors typically provide accurate 
measurements at high sampling rates. The interpolation 
is done using a Bayesian non-parametric approach where 
prior model knowledge is incorporated in a Gaussian 
process (GP) prior. GPs (Rasmussen and Williams 2006) 
are powerful tools for Bayesian non-parametric inference 
and learning, and they provide a framework for fusing 
first-principles prior knowledge with quantities of noisy 
data. 

The contributions of this paper are three-fold. Eirst, 
we model the ambient magnetic field using a Gaussian 
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process prior in which we incorporate physical knowledge 
about the magnetic field. This extends the work by 
Wahlstrom et al. (2013) by presenting an approach where 
the GP prior is a latent (unobservable) magnetic potential 
function. Second, we use a computationally efficient GP 
implementation that allows us to use the large amounts 
of data provided by the magnetometer. To circumvent the 
well-known computational challenges with GPs (see, e.g., 
Rasmussen and Williams 2006), we rewrite the model 
in terms of a Hilbert space representation introduced by 
Solin and Sarkka (2014). We extend the approach to allow 
for modeling of the bias caused by the Earth magnetic 
field. Third, we use this method in combination with the 
sequential approach introduced in Sarkka et al. (2013). 
This allows for online updating of the magnetic field 
estimate. It also opens up the possibility to focus on the 
spatio-temporal problem in which the magnetic field can 
change over time, for instance due to furniture being 
moved around. The proposed method is experimentally 
tested on a number of data sets, showing that accurate 
estimates of the magnetic field can be obtained. 

Our interest in interpolation of the magnetic field 
is for future use in indoor positioning and navigation 
applications. In these applications, measurements that are 
accurate on a short time-scale—but drift on longer time- 
horizons—are typically combined with absolute position 
information. Examples of this use data from wheel 
encoders or from inertial sensors together with ultra- 
wideband, Wi-Ei, or optical measurement equipment 
such as cameras (see, e.g.. Woodman 2010; Hoi 2011). 
The downside of these sources of absolute position is 
that they typically rely on additional infrastructure or 
require certain conditions to be fulfilled such as line- 
of-sight. The advantage of using the magnetic field for 
positioning is that it can be measured by a small device, 
without additional infrastructure and without line-of-sight 
requirements. Eurthermore, magnetometers are nowadays 
present in (almost) any inertial measurement unit (IMU) 
or smartphone. A requirement for localization using the 
ambient magnetic field as a source of position information 
is that accurate maps of the magnetic field can be 
constructed within reasonable computational complexity. 
This can also be regarded as a step towards simultaneous 
localization and mapping (SLAM) using magnetic fields 
in which localization is done while building the map (see. 


e.g., Leonard and Durrant-Whyte 1991; Durrant-Whyte 
and Bailey 2006). 

This paper is structured as follows. The next section 
covers a survey of existing work, which also provides 
additional motivation for the approach. Section 3 gives the 
necessary and sufficient background of magnetic fields, 
while a more detailed walk-through is provided in the 
appendices. The Gaussian process regression model is 
constructed in Section 4, which is then extended to 
explicit algorithms for batch and sequential estimation 
in the next section. Section 6 covers several different 
experiments. The experiments and some additional 
comments regarding the methodology are discussed at the 
end of the paper. 

2 Related work 

Spatial properties of the magnetic field have been of 
interest in a large variety of research domains. Eor 
instance, the magnetic field has been extensively studied 
in geology (see, e.g., Nabighian et al. 2005) but also in 
magnetospheric physics, geophysics, and astrophysics. In 
all of these domains, interpolation of the magnetic field is 
of interest (see, e.g., Guillen et al. 2008; Calcagno et al. 
2008; Mackay et al. 2006; Bhattacharyya 1969; Springel 
2010, for examples of magnetic field interpolation in the 
respective areas). 

In recent years interest has emerged in using the 
magnetic field as a source of position information for 
indoor positioning (Haverinen and Kemppainen 2009). 
Eeasibility studies have been conducted, focusing both 
on the time-varying nature of the magnetic field and on 
the amount of spatial variation in the magnetic field. Li 
et al. (2012) report experiments showing that the magnetic 
field in a building typically shows large spatial variations 
and small time variations. This is also supported by the 
experimental study reported by Angermann et al. (2012) 
in which significant anomalies of the ambient magnetic 
field are reported. These experiments give confidence 
that the magnetic field provides sufficient information 
for localization purposes. However, Li et al. (2012) also 
report significant temporal changes in the magnetic field 
in the vicinity of mobile magnetic structures, in their case 
an elevator. 

A number of approaches have been reported on 
building a map of the ambient magnetic field for 
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indoor localization purposes. Le Grand and Thrun 
(2012) propose a method to build a map of the 
magnetic field by collecting magnetometer data in a 
grid and linearly interpolating between these points. 
This map is subsequently used for localization with a 
particle filter combining magnetometer and accelerometer 
measurements from a smartphone. Robertson et al. (2013) 
present a SLAM approach for pedestrian localization 
using a foot-mounted IMU. They use the magnetic 
field intensity which they model using spatial binning. 
Frassl et al. (2013) discuss the possibility of using more 
components of the magnetic field (for instance the full 
three-dimensional measurement vector) in the SLAM 
approach instead. Vallivaara et al. (2010, 2011) present 
a SLAM approach for robot localization. They model 
the ambient magnetic field using a squared exponential 
GP prior for each of the magnetic field components. 
Wahlstrom et al. (2013) incorporate additional physical 
knowledge by making use of Maxwell’s equations 
resulting in the use of curl- and divergence-free GP priors 
instead. 

As can be concluded, there exists a wide range 
of existing literature when it comes to modeling the 
ambient magnetic field. The amount of information that 
is used differs between the approaches. For instance, 
some approaches use full three-dimensional magnetic 
field vectors while others only use a one-dimensional 
magnetic field intensity. Furthermore, the amount of 
physical information that is included differs. In this 
paper, we build on the approach by Wahlstrom et al. 
(2013) and use the full three-dimensional magnetometer 
measurements. We include physical knowledge in terms 
of the magnetic field potential. 

As discussed above, GPs have frequently been used 
in modeling and interpolation of magnetic fields. GP 
regression has also successfully been applied to a wide 
range of applications (see, e.g., O’Callaghan and Ramos 
2012; Smith et al. 2011). Furthermore, it has previously 
been used for SLAM (see, e.g., Tong et al. 2013; Ferris 
et al. 2006; Barkby et al. 2012; Barfoot et al. 2014; 
Anderson et al. 2015). One of the challenges in using GPs 
is the computational complexity (see, e.g., Rasmussen and 
Williams 2006), which scales cubicly with the number of 
training data points. Considering the high sampling rate 
of the magnetometer and the fact that each observation 
contains three values, a large number of measurements 
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Earth magnetic field Ferromagnetic object Distorted field 

Figure 1. A ferromagnetic object deflects the Earth’s 
magnetic field and introduces distortions in the field. 



is typically available for mapping. Because of these 
computational challenges, the data in Wahlstrom et al. 
(2013) was downsampled. 

Attempts to speed up GP inference have spawned 
a wide range of methods which aimed at bringing 
GP regression to data-intensive application fields. These 
methods (see Quinonero-Candela and Rasmussen 2005, 
for a review) typically build upon reducing the rank 
of the Gram (covariance) matrix and using the matrix 
inversion lemma to speed up matrix inversion. For 
stationary covariance functions, the spectral Monte Carlo 
approximation by Lazaro-Gredilla et al. (2010), or the 
Laplace operator eigenbasis based method introduced by 
Solin and Sarkka (2014) can be employed. For uniformly 
spaced observations, fast Fourier transforms can provide 
computational benefits (Paciorek 2007; Fritz et al. 2009). 
As will be shown later on in this paper, the Laplace 
operator approach by Solin and Sarkka (2014) falls 
natural to modeling of the magnetic field in terms of a 
magnetic field potential. 

All approaches on mapping of magnetic fields 
discussed above assume that the magnetic field is constant 
over time. However, as shown by Li et al. (2012) 
significant temporal changes in the magnetic field occur 
in the vicinity of mobile magnetic structures. For GP 
models evolving in time, spatio-temporal GP models 
can be solved efficiently using Kalman filtering methods 
(Hartikainen and Sarkka 2010; Osborne 2010; Sarkka 
et al. 2013; Huber 2014). In this paper, we will take the 
approach of Sarkka et al. (2013) to compose a spatio- 
temporal GP prior for the model and solve the inference 
problem by a sequential Kalman filtering setup. This 
allows for online estimation of the magnetic field estimate 
and can be used to allow for time variations in the 
magnetic field. 
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3 The ambient magnetic fieid 

On a macroscopic scale, magnetic fields are vector fields, 
meaning that at any given location, they have a direction 
and strength (magnitude). These properties are familiar 
from everyday life: the force created by permanent 
magnets attracting and repelling ferromagnetic materials 
is used in various utensils, and the compass aligning itself 
with the direction of the Earth’s magnetic field has proved 
invaluable for mankind during the past centuries. The 
Earth’s magnetic field sets a background for the ambient 
magnetic field, but deviations caused by the bedrock and 
anomalies induced by man-built structures defiect the 
Earth’s magnetic field. This makes the magnetic field vary 
from point to point, see Eigure 1. 

We describe the magnetic field with a function H(x), 
where H : ^ Eor each point in space x there will 

be an associated magnetic field H(x). Such a vector field 
can be visualized with field lines, where points in space 
are associated with arrows. The principles under which 
the magnetic field is affected by structures of buildings 
are well known and governed by the very basic laws of 
physics, see Appendix A for more details. 

By further following Appendix A, the magnetic field H 
is shown to be curl-free, meaning that 

VxH = 0. (1) 

Eor this relation to hold we assume that there is no 
free current (current in wires for example) in the region 
of interest. Such an assumption is valid in most indoor 
environments where the major source for variations in the 
ambient field is caused by metallic structures rather than 
free currents in wires. 

One property of curl-free vector fields is that the line 
integral along a path P only depend on its starting point 
A and end point B, and not on the route taken 

J H(x) ■ dx = (p(A) - (2) 

where ip scalar function ^ M. This can be 

rewritten by interpreting cp sls sl scalar potential 

H = -Wcp. (3) 

In Eigure 2 illustrates the curl-free property and the scalar 
potential. Domain ft is curl-free and has an associated 



Figure 2. Illustration of a vector field with non-zero curl. The 
vortex point makes it non-curl-free as the vector field curls 
around it. However, the subset ft excludes the vortex point 
and the vector field is curl-free in this region. To this region a 
scalar potential p can be associated, here illustrated with 
shading. 

scalar potential, while the entire domain is not curl-free 
due to the vortex point. Note that in a non-curl-free 
vector field no such scalar potential exists since a line 
integral around the swirl is non-zero. Eor the magnetic 
field H, the swirl corresponds to a wire of free current 
pointing perpendicular to the plane, which we assume is 
not included in the region of interest. 

The relation (3) is the key equation that we will exploit 
in our probabilistic model of the ambient magnetic field. 
We will choose to model the scalar potential p instead 
of the magnetic field H directly. This implicitly imposes 
the constraints on the magnetic field that the physics 
is providing. This model will be explained in the next 
section. 

4 Modeling the magnetic field using 
Gaussian process priors 

In this section, we introduce our approach to modeling 
and interpolation of the ambient magnetic field. We 
use a Bayesian non-parametric model in which we use 
knowledge about the physical properties of the magnetic 
field as prior information. We tackle the problem of 
interpolating the magnetic field using Gaussian process 
(GP) regression. In Section 4.1 we first give a brief 
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background on GPs. After this, we introduce the problem 
of modeling the magnetic field in Section 4.2. A 
commonly used GP model of the magnetic field is 
introduced in Section 4.3. In Section 4.4, we subsequently 
introduce our proposed method for modeling the magnetic 
field in which we encode the physical properties from 
magnetostatics that were presented in the previous 
section. 

4. / Gaussian process regression 

In GP regression (Rasmussen and Williams 2006) the 
model functions /(x) are assumed to be realizations from 
a Gaussian random process prior with a given covariance 
function /^(x, x'). Learning amounts to computing the 
posterior process at some test inputs x>^ given a 
set of noisy measurements , ^ 2 , • • •, observed at 
xi, X 2 ,..., Xn, respectively. This model is often written 
in the form 

/(x) ~ aP(0,K(x,x')), 

Vi = /(Xi) +£*, 

where the observations yi are corrupted by Gaussian 
noise Si ^ N(0, cr^oise)’ for i = 1, 2,..., n. Because both 
the prior and the measurement model are Gaussian, 
the posterior process will also be Gaussian. Hence, the 
learning problem amounts to computing the conditional 
means and covariances of the process evaluated at the test 
inputs. 

Prediction of yet unseen process outputs at an input 
location x* amounts to the following in GP regression: 
p(/(x*) I V) = N(/(x*) I E[/(x*)], V[/(x*)]). The 
conditional mean and variance can be computed in 
closed-form as (see Rasmussen and Williams 2006, here 
the conditioning is left out in the notation for brevity) 

E[/(x,)]=kI(K + aL,I„)-V, 

V[/(x*)] = k(x*,x*) -kJ(K + c 72 „i 3 gI„)-ik*, 

where 'Kij = k* is an n-dimensional vector 

with the ith entry being K(x^,Xi), and y is a vector of 
the n observations. Furthermore, due to Gaussianity, the 
marginal likelihood (evidence) of the covariance function 
and noise parameters can also easily be computed, 
allowing for Bayesian inference of the parameters as well 
(Rasmussen and Williams 2006). 


The choice of a specific covariance function encodes 
the a priori knowledge about the underlying process. 
One of the most commonly used covariance functions, 
which will also frequently be used in the next sections, 
is the stationary and isotropic squared exponential (also 
known as exponentiated quadratic, radial basis function, 
or Gaussian). Following the standard notation from 
Rasmussen and Williams (2006) it is parametrized as 

kse(x, x') = (t|e exp ^ ^ , (6) 

where the hyperparameters alg and ^se represent the 
magnitude scale and the characteristic length-scale, 
respectively. These can be learned from data, for instance 
by maximizing the marginal likelihood. 

4.2 Interpolation of magnetic fields 

In this work we tackle the problem of interpolating the 
magnetic field to spatial locations from where we do not 
have any measurements. In other words, we will tackle 
the problem of predicting the latent (unobservable noise- 
free) magnetic field f(x>,) (such that f : ^ M^) at 

some arbitrary location x* given a set of noise-corrupted 
measurements = {(x^, of the magnetic field. 

Here, the measurements y correspond to the H-field 
corrupted by i.i.d. Gaussian noise. 

Two important things need to be noted with regard to 
the interpolation of the magnetic field. First, note that 
the measurements of the magnetic field are vector-valued 
(contrary to the scalar observations in (4)). This raises 
the question of how to deal with the different magnetic 
field components. They can either be treated separately 
as will be done in Section 4.3, or a relation between the 
different components can be assumed, as is the case in 
the method we propose in Section 4.4. Secondly, note 
that the function describing the magnetic field is not zero- 
mean, contrary to the GP model in (4). Instead, its mean 
lies around a local Earth magnetic field. This depends 
on the location on the Earth but can also deviate from 
the Earth’s magnetic field in indoor environments due to 
magnetic material in the structure of the building. The 
unknown mean can be modeled as an additional part of the 
covariance function /^(x, x') (Rasmussen and Williams 
2006). 
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A B C D ED E 

(a) The route (b) The data and the GP prediction for D-E 

Figure 3. A simulated example of the interpolation problem, (a) Training data has been collected along the route A-D, but 
the magnetic field between D-E is unknown, (b) The noisy observations of the magnetic field between A-D, and GP 
predictions with 95% credibility intervals. Both the independent GP modeling approach (with shared hyperparameters) and 
the scalar potential based curl-free GP approach are visualized. The simulated ground truth is shown by the solid lines. 


An illustration of GP regression for magnetic fields is 
provided in Figure 3, where noisy readings of a magnetic 
field have been collected along route A-D (comprising 
V), and the magnetic field along route D-E (comprising 
the prediction locations x*) needs to be inferred from 
the measurements. Each component varies around a local 
magnetic field due to magnetic material in the vicinity 
of the sensor. Different interpolation techniques can be 
used based on different prior knowledge that can be 
incorporated in the GP. Two different interpolation results 
are shown: one based on independent modeling of each 
vector field components (with shared hyperparameters) 
and another based on associating the GP prior with the 
scalar potential of the (curl-free) vector field. These are 
based on the models we will introduce in the coming two 
sections. 

4.3 Separate modeling of the magnetic field 
components 

The most straightforward approach to GP modeling of 
vector-valued quantities is to model each of the field 
components as an independent GP. This approach has 
been widely applied in existing robotics literature (see, 
e.g ., Vallivaara et al. 2010, 2011; Kemppainen et al. 2011; 
Jung et al. 2015; Viseras Ruiz and Olariu 2015). For each 


of the three magnetic field components d G {1, 2, 3}, this 
model can be written as 

/d(x) -- t/P(0, Kconst.(x,x') + Kse(x,x')), 

yd,i = fdi^i) + ei.d, 

where the observations i/d^, i = 1, 2,..., n, are corrupted 
by independent Gaussian noise with variance cr^oise- 
non-zero mean of the magnetic field is handled by a 
constant covariance function 

^const. (^5^ ) = ^const. 5 (^) 

where cr^onst. ^ magnitude scale hyperparameter. The 
small-scale variation in the field is modeled by a squared 
exponential covariance function (6). Hence, the model (7) 
encodes the knowledge that the realizations are expected 
to be smooth functions in space with a constant shift from 
zero mean. 

The model has four hyperparameters: two magnitude 
scale parameters (cr^onst. ^se)’ ^ length-scale 

parameter (^se), and a noise scale parameter (cr^oise)- 
Assuming that the components are completely separate, 
each component has four hyperparameters to learn. The 
resulting model is fiexible, as it does not encode any 
relation between the vector field components. In practice, 
this might lead to problems in hyperparameter estimation. 
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with parameter estimates converging to local optima and 
magnetic field components behaving very differently with 
respect to each other. Therefore, the hyperparameters are 
often fixed to reasonable values—instead of learned from 
data. 

A more sensible approach for separate, but not 
completely independent, modeling of the magnetic field 
measurements, models them as realizations of three 
independent GP priors with joint learning of the shared 
hyperparameters (see also Kemppainen et al. 2011). Note 
that for this model, the covariance in the GP posterior is 
independent of the outputs y and only depends on the 
input locations x (which are shared for all components 
in y). Hence, calculating the marginal likelihood only 
requires inverting a matrix of size n (not 3n). For this 
model, the expression for evaluating the log marginal 
likelihood function for hyperparameter optimization can 
be written as 

m = - logp(y \e,V) = ^- log IKe + aL, I„| 

1 Sr? 

+ 2+ ^noise In)” V^] + y log(27r), (9) 

where y G and Even though this 

formulation where the matrix inversion is only of size 
n is not a new trick, this formulation of joint learning 
of hyperparameters is apparently not commonly used in 
robotics applications. 

Figure 3b shows the results of predicting the magnetic 
field behavior along the route D-E for the GP prior 
modeling the three magnetic field components separately 
but with joint learning of the shared hyperparameters. 
The colored patches show the 95% credibility intervals 
for the prediction with the mean estimate visualized by 
the white line. The simulated ground-truth (solid colored 
line) falls within the shown interval, and the model 
captures the general shape of the magnetic field variation 
along the path. The strengths of this model are that it is 
flexible and that the assumptions are conservative. The 
weaknesses on the other hand are evident: The model 
does not incorporate physical knowledge of the magnetic 
field characteristics. In the next section we will instead 
explore this knowledge by modeling the magnetic field as 
derivative measurements of a scalar potential. 


4.4 Modeling the magnetic field as the 
gradient of a scalar potential 

Following our choices in Section 3, we assume that the 
magnetic field H can be written as the gradient of a 
scalar potential (p(x) according to (3). Here, (/? : ^ M 

and X G is the spatial coordinate. We assume (/?(x) 
to be a realization of a GP prior and the magnetic field 
measurements y^ G to be its gradients corrupted by 
Gaussian noise. This leads to the following model 

(p(x) -- gv{0, /^lin.(x, x') + /^se(x, x')), 

. M ( 10 ) 

yi = -v¥^(x)lx=x,+ 

where ^ N( 0 , I 3 ), for each observation i = 

1, 2,... ,?i. Recall that nabla (V) is a linear operator, 
and that Gaussianity is preserved under linear operations. 
Thus, we can still derive a closed-form solution for 
this GP model. For simplifying the notation in the next 
sections, we introduce the notation f (x) for the gradient 
field evaluated at x. 

The squared exponential covariance function ( 6 ) in (10) 
allows us to model the magnetic field anomalies induced 
by small-scale variations and building structures. The 
local Earth’s magnetic field contributes linearly to the 
scalar potential as 

Klin.(x,x') = (11) 

where is the magnitude scale hyperparameter. The 
derivative of this linear component contributes to the 
model in a similar way as the constant covariance function 
( 8 ) in the previous section. 

Our model now has four hyperparameters: two 
magnitude scale parameters and cr|g), a length-scale 
parameter (^se), and a noise scale parameter (cr^oise)- We 
learn these parameters from the data, by maximizing the 
marginal likelihood. 

The second set of predictions for the route D-E in 
Figure 3b shows the interpolation outcome from using 
the above model. In comparison to the independent GP 
model, the scalar potential based GP prior provides 
additional information to the model by tying the vector 
field components to each other. This improves the 
estimates in terms of accuracy and makes the 95% 
credibility interval more narrow. 
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Incorporating physical knowledge about the magnetic 
field into a GP prior is not new. In fact, the assumptions 
in our model ( 10 ) are similar to the assumptions made by 
Wahlstrom et al. (2013). Their assumptions are derived 
from Maxwell’s equations and model the magnetic field 
as curl-free using a multi-output GP 

f (x) ~ gv{0, I3 + Kcuri(x, x')), 
y* =f(xi) + £i, 

such that f : ^ m 3 and Si ~ N( 0 ,(T 2 ^i 3 j 3 ). The 

curl-free covariance function Kcuri : x ^ is 

given by 


Kcuri (x 7 X ) — 





x') (x-xyi 


exp 


l|x-x-|n 

) ' 
(13) 


Note that the parametrization used by Wahlstrom et al. 
(2013) is slightly different in terms of the definitions 
of the hyperparameters. The curl-free kernel (13) 
was discussed in, for example, Fuselier (2007) and 
Baldassarre et al. (2010). The topic is also covered by 
Alvarez et al. ( 2012 ). 

As shown in Appendix B, our proposed model (10) 
and the model (12) are actually equivalent. In this 
work, however, the model formulation through the scalar 
potential is crucial, as it will enable us to easily extend 
the model to an approximate form for efficient GP 
inference in the next section. We will, however, report the 
magnitude scale hyperparameter values as given through 
the relation to the curl-free covariance function: 

This ensures the quantity to have more meaningful units 

(/iT). 


approach more or less useless in practice, when the 
number of observations becomes large (say n > 1000). 

This is a fundamental restriction associated with the 
naive formulation of GP models involving the inversion 
of the covariance matrix. Using special structure of 
the problem and/or approximative methods, this high 
computational cost can often be circumvented. In this 
section we present an approach which both uses the 
special differential operator structure and projects the 
model on a set of basis functions characteristic to the 
covariance function. We first present the method for 
spatial batch estimation, and then extend it to a temporal 
dimension as well. 

Existing GP methods for mapping and interpolation 
of the magnetic field have been considering only batch 
estimation, where the data is first acquired and then 
processed as a batch. In this section, we aim to extend this 
to an online method, enabling the GP regression estimate 
of the magnetic field to be updated when new data is 
acquired. We denote such a dataset as Vn = {(xi,yi) | 
i = 1, 2,..., n}, and thus Vi denotes all the data that has 
been observed up to time instance ti . 

Considering time as part of the data stream enables us 
to think of three distinctive setups for estimation of the 
magnetic field: 

• Batch estimation of the magnetic field, where the 
data is first acquired and then the field is estimated at 
once. 

• Sequential updating of the field estimate, where we 
assume all the measurements to be of the same static 
magnetic field. 

• Spatio-temporal estimation of the time-dependent 
magnetic field, where we assume the field to change 
over time. 


5 Efficient GP modeling of the magnetic 
field 

Gaussian processes are convenient tools for assigning 
fiexible priors to data—as we saw in the previous 
section. However, the main problem with the model in 
the previous section is its high computational cost. The 
approach scales as 0{n^) (recall that each observation 
is three-dimensional, meaning 3n becomes large very 
quickly). This computational complexity renders the 


In the next sections we will present how these 
scenarios can be combined with the scalar potential 
based GP scheme without requiring to repeat the batch 
computations after each sample. 

5.1 Reduced-rank GP modeling 

A recent paper by Solin and Sarkka (2014) presents an 
approach that is based on a series expansion of stationary 
covariance functions. The approximation is based on the 
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following truncated series: 

m 

k(x,x') « 0j(x'), (14) 

i=i 

where 5'(-) is the spectral density of the covariance 
function /^(•, the jth eigenfunction of the 

negative Laplace operator and is the corresponding 
eigenvalue. The efficiency of this approach is based on 
two properties: (i) the eigenfunctions are independent 
of the hyperparameters of the covariance function, 
and (ii) for many domains the eigenfunctions and 
eigenvalues can be solved beforehand in closed-form. 
Truncating this expansion at degree m n allows the 
GP regression problem to be solved with a 0{nrin?) and 
the hyperparameters to be learned with a 0{m^) time 
complexity. The memory requirements scale as 0{nm). 

However, what makes this approach even better suited 
for this problem, is that the approximation is based on 
the eigendecomposition of the Laplace operator. This 
eigenbasis falls natural to the problem formulation in 
the previous section, where the latent potential field 
is observed through gradients. As for the choice of 
covariance structure, the squared exponential covariance 
function is stationary, so there is no problem applying this 
approach for that part of the above problem formulation. 
The linear covariance function is not stationary, but as we 
will see, that is not a problem in this case. 

Our interest lies in modeling the magnetic field in 
compact subsets of allowing us to restrict our 
interest to domains Q comprising three-dimensional 
cuboids (rectangular boxes) such that xG [—Li^Li] x 
[— L 2 , 1 / 2 ] X [— 1 / 3 , 1 / 3 ] C (recall that a stationary 
covariance function is translation invariant). In this 
domain, we can solve the eigendecomposition of the 
Laplace operator subject to Dirichlet boundary conditions 

f-V2(/)j(x) = A20 ^(x), xen, 

I 0i(x) = 0, X G dfl. 

The choice of the domain and boundary conditions is 
arbitrary, but for regression problems with a stationary 
covariance function the model reverts back to the prior 
outside the region of observed data, so the Dirichlet 
boundary condition does not restrict the modeling if ft is 
chosen suitably. 


This particular choice of domain and boundary 
conditions yield the following analytic expression for the 
basis functions: 



where the matrix n G consists of an index 

set of permutations of integers (i.e., 

{(1,1,1), (1,1,2),..., (1, 2,1),..., (2,1,1),...}). The 
basis functions are independent of the hyperparameters, 
and thus only need to be evaluated once. Adding the 
linear covariance function to the GP prior corresponds 
to inserting a set of three linear basis functions into the 
model. 

The computational benefits come from the approximate 
eigendecomposition of the Gram (covariance) matrix, 
=/^(xi,Xj) (see Solin and Sarkka 2014, for 
derivations and discussion). It can now be written out 
in terms of the basis functions and spectral densities: 
K ^ ^A^^. The basis functions which span the solution 
are collected in the matrix ^ with the 

following rows 

= (x7,0l(Xi),02(Xi),---,0m(Xi)) , (18) 

for i = l,2,...,n. Accordingly, we define the cor¬ 
responding measurement model matrix projecting the 
derivative observations onto the basis functions. Analo¬ 
gously, we define the matrix G ]^3nx(3+m) 
following block-row matrix: 

= (Vx7,V0i(Xi),V(/)2(Xi),...,V(/>m(Xi)) , 

(19) 

for i = 1, 2,..., n. Similarly we define and 
as vectors evaluated at the prediction input location 
X* defined analogously to Equations (18) and (19), 
respectively. The diagonal matrix A is defined by 

A = diag((T|„ , cr,i , ai„ , 

5'se(A 2), . . . , 5'sE(Am))- (20) 

For three-dimensional inputs, the spectral density 
function of the squared exponential covariance function 
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(6) is given by 

S'se(w) = (tIe {2'K£l-^f/‘^ , ( 21 ) 

where the hyperparameters and ^se characterize the 
spectrum. 

The Laplace operator eigenbasis approximation 
method can be combined with the independent 
GP approach, the independent GPs with shared 
hyperparameters, and the scalar potential GP approach. 
In the next sections the presentation will be specific to 
the scalar potential model, but a similar setup can be 
constructed for the other methods as well. 

5.2 Batch estimation 

We first tackle the batch estimation problem which 
provides the approximative solution to the GP regression 
problem in Equation (5) for the scalar potential GP. 

Following the derivations of Solin and Sarkka (2014), 
predictions for interpolation and extrapolation of the 
magnetic field at yet unseen input locations x* are given 
by: 

IE[f(x*)] 

« V#*([V#]Tv# + <TL,,A-i)-i[V#]Tvec(y), 
V[f(x,)] 

^ <ise V#*([V#]TV# + 

( 22 ) 

where vec(-) is the vectorization operator which converts 
a matrix to a column vector by stacking its columns on 
top of each other, such that the 3 x n matrix is converted 
into a vector of size 3n. The basis functions and 
need to be evaluated by Equation (19), and A by (20). 

For this model, the expression for evaluating the 
log marginal likelihood function for hyperparameter 
optimization can be written as 

C{e) = ^log\Ke + aiM 

1 3n 

+ 2 vec(y)’^(K0 + l3„)“^vec(y) + y log(27r), 

( 23 ) 


Algorithm 1. Algorithm for batch estimation of the 
scalar potential GP magnetic field with the reduced-rank 
approach. 

Input: V = {(xj, yi)}"=i, x*, fl, m. 

Output: E[f(x*)],V[f(x*)]. 

1: Use Eq. (19) to evaluate the basis functions 
from x^s and U. 

2: Use Eq. (23) to optimize hyperparameters 

3: Use Eq. (19) to evaluate the basis functions 
from x^s and U. 

4: Solve the GP regression problem by Eq. (22). 


where the quantities can be approximated by 
log IKe + l 3 „| « (3n - m) log 

m 

+ + log l^noiseAe ' + [V#]Tv$|, (24) 

i=i 

vec(y)’'’(K0 + 

i^noise Isn) ^vec(y) 

~ [vec(y)'^vec(y) - vec(y)'^V#(cr2„iseA0^ 

^ noise 

+ [V#]'rV$)-MV$]'rvec(y)], (25) 

where the only remaining dependency on the covariance 
function hyperparameters are in the diagonal matrix A 
defined through the spectral density in Equation (20). In a 
software implementation, Cholesky decompositions can 
be employed for numerical stability in the calculation 
of determinants and matrix inverses. For optimizing 
the hyperparameters, gradient based optimizers can be 
employed (see Solin and Sarkka 2014, for details on 
deriving the partial derivatives). 

Algorithm 1 describes the step-by-step workfiow for 
applying these equations in practice. The inputs for the 
method are the data V (spatial points and the magnetic 
field readings), the test points x* to predict at, the 
domain boundaries U, and the approximation degree 
parameter m (controlling the accuracy of the Hilbert 
space approximation, see Eq. (14)). The algorithm returns 
the marginal mean and variance of the predicted magnetic 
field at X*. The scalar potential could be returned instead 
by using Equation (18) in step three. 
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5.3 Sequential estimation 

In robotics, considering an online solution for estimating 
the magnetic field estimates falls natural. The following 
formulation provides the same (within numerical preci¬ 
sion) solution as the batch estimation solution ( 22 ) in the 
previous section. The inference scheme in the previous 
section is in practice the solution of a linear Gaussian 
estimation problem. Sequential solutions for this type of 
problems have been extensively studied, and this mathe¬ 
matical formulation is widely known as the Kalman filter. 
The connection between Kalman filtering and Gaussian 
process regression has recently been studied, for example, 
by Osborne (2010); Sarkka et al. (2013); Huber (2014). 

Reformulation of a batch problem to a sequential 
algorithm is discussed (with examples) in the book Sarkka 
(2013). Following this formulation (and notation to large 
extent) we may write the following recursion: Initialize 
/L^o = 0 and 5]o = (from the GP prior). For each new 
observation i = 1 , 2 ,..., n update the estimate according 
to 

Si = V^i5]i_i[V^i]"'" -f cr^oise 

= Hi-i + Kj(yj - 

Si = Si_i - KiSiKj. 

This means that for a test input location x>^ we get 
predictions for the mean and the variance of the magnetic 
field which are given by 

E[f(x*) I Vi] « V#* ^li, 

V[f(x*) I Vi] « V#* Si 

and conditional on the data observed up to observation 
i. Writing out the conditioning on V was stripped in the 
earlier sections for brevity. 

Here we do not consider optmization of the hyper¬ 
parameters 6. The marginal likelihood can be evaluated 
through the recursion, but in an online setting we suggest 
optimizing the hyperparameters with some initial batch 
early in the data collection and then re-optimizing them 
later on if necessary. 

Algorithm 2 presents the scheme of how to apply the 
equations in practice. The inputs are virtually the same 
as for the batch algorithm, but now the hyperparameters 


Algorithm 2. Algorithm for sequential modeling of the 
scalar potential GP magnetic field estimate. Alternative 
(a) corresponds to the sequential model, and (b) to the 
spatio-temporal modeling approach. 

Input: Vn, X*, ft, m, 6. 

Output: Bar. 

1: Initialize fio = 0 and 5]o = from Eq. (20). 

2: for i = 1, 2,..., n do 

3: Evaluate by Eq. (19) from x^. 

4: (a) Perform an update by Eq. (26). 

(b) Perform an update by Eqs. (31-32). 

5: Evaluate the current prediction at x* by Eq. (27). 

6: end for 


6 are considered known. The current estimate can be 
returned after each iteration loop. 

5.4 Spatio-temporal modeling 

The sequential model allows for extending the modeling 
to also track dynamic changes in the magnetic field 
without virtually any additional computational burden. 
Let the data Vn = {(ti, x^, yi)}^=i now also comprise a 
temporal variable t which indicates the time when each 
observation was acquired. 

We present the following spatio-temporal model for 
tracking changes in the ambient magnetic field. The 
spatio-temporal GP prior assigned to the scalar potential 
(/?(x, t), depending on both location and time, is defined 
as follows: 

gvio, Kii„,(x,x') +KsE(x,x')Kexp(i,^'))) 

(28) 

where the additional covariance function /^exp(f,fO 
defines the prior assumptions of the temporal behavior. 
This covariance function is defined through 

«^exp = exp( - J, (29) 

where ^time is a hyperparameter controlling the length- 
scale of the temporal effects. In the temporal domain, 
this model is also known as the Ornstein-Uhlenbeck 
process (see, e.g., Rasmussen and Williams 2006). The 
assumption encoded into it is that the phenomenon is 
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continuous but not necessarily differentiable. Therefore it 
provides a very flexible means of modeling the changing 
ambient magnetic held. Also note that in Equation (28) 
the temporal effects are only associated with the anomaly 
component, the bias being tracked as a static component. 

Following the derivations in Hartikainen and Sarkka 
( 2010 ), we can write down the dynamic state space model 
associated with the time evolution of the spatio-temporal 
GP prior model (28) 

Aj = blkdiag(l3, exp(-Atj/ 4 n,e)), 

Qi = blkdiag(03,Im[l - exp(-2Ati/4nie)]), 

where Ati = — ti is the time difference between two 

consecutive samples and O 3 denotes a 3 x 3 zero matrix. 

For the time update (Kalman prediction step) we may 
thus write 


+ Qi-l, 


(31) 


and the modifled measurement update (Kalman update 
step) 


p,i = ij.i + Ki(yi - 
Si = Si - KiSiKj. 


(32) 


Algorithm 2 features the workflow of applying the 
method (option (b)). In practice the only additional input 
is including the temporal length-scale in 6. 


6 Experiments 

The experiments in this paper are split into four parts. 
We first demonstrate the feasibility of the scalar potential 
approach with simulated data, where comparison to 
known ground truth is possible. After this we present 
a small-scale proof-of-concept demonstration of the 
approach. 

The final two examples are closer to real-world use 
cases. The third experiment is concerned with mapping 
the magnetic held in a building using a handheld 
smartphone. The final experiment uses an inexpensive 
mobile robot for online mapping and real-time tracking 
of the changing magnetic held. 


6.1 Simulated experiment 

As a first part of our experimental validation of the 
model, we present a simulation study. This will be used 
to illustrate our method and to quantify its performance 
through Monte Carlo simulations. In the simulations, 
we assume that the magnetic held measurements can 
indeed be modeled using a scalar potential as argued in 
Section 3. Hence, we simulate the magnetometer data 
from the GP that models the magnetometer measurements 
as gradients of a scalar potential held as discussed in 
Section 4.4. To reduce the computational complexity, 
we simulate the data using the computationally efficient 
approach described in Section 5.2 with a large number of 
basis functions (m = 4096). This reduces computational 
complexity and is, as shown by Solin and Sarkka (2014), 
a good approximation of the true model. 

The magnetometer is assumed to move in a three- 
dimensional volume with x^y^z G [—0.4, 0.4] such that 
X = {x^y^z). The training data used for training the GP 
is randomly uniformly distributed over this volume. A 
validation data set is used to assess the predictive power of 
the trained GP. This validation data is a three-dimensional 
meshgrid over the same volume and consists of nyai. = 
9261 positions Xyai. with a true magnetic held ftrue(xvai.)- 
The magnetic held is predicted at these points using 
the trained GP, leading to ftrain.(xvai.), after which the 
quality of the GP solution can be assessed in terms of 
the root mean square error (RMSE). We set the domain 
Q in the GP model to Li = 1/2 = Ls = 0.5 and simulate 
using the following hyperparameters cr^onst. =0.3, cr|g = 
1 ,^SE = 0 . 1 , and 4.4 for 

more details on the notation). 

We compare our proposed method with two other 
approaches. The first is the approach from Vallivaara 
et al. ( 2010 , 2011 ), where the magnetic held is modeled 
using independent GPs for all three components. Each GP 
consists of a constant and a squared exponential kernel 
and has its own hyperparameters. The second approach, 
considered in Kemppainen et al. (2011), models the 
magnetic held similarly but with shared hyperparameters. 
For details, see also Section 4.3. 

In a first set of Monte Carlo simulations, we analyze the 
performance of the three different approaches depending 
on the number of (randomly distributed) training data 
points, that is in terms of the sparseness of the 
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Figure 4. The average RMSE (with standard deviations) from 30 Monte Carlo simulations as a function of the number of 
simulated data points used for training the GR The constantly lower error for the Hilbert space scalar potential GP in 
comparison to the separate GP models and shared hyperparameter GP model is explained by the additional prior physical 
knowledge encoded into the model. In (b), the shared hyperparameter GP shows a slight advantage over the fully 
independent models. 


magnetometer data and the amount of interpolation that 
is needed for prediction. For all three approaches we 
use a large number of basis functions (m = 4096). They 
are hence expected to approach the performance of the 
full GP solution. To exclude problems with local minima 
in the hyperparameter optimization—which will be the 
topic of a second set of simulations—the hyperparameter 
optimization is started in the values used for simulating 
the data. The results from 30 Monte Carlo simulations 
are shown in Figure 4a. Naturally, the more data is used 
for training the GP, the smaller the RMSE becomes. 
The GP that models the magnetic field measurements as 
gradients of a scalar potential field, outperforms the other 
two approaches, independent of the amount of training 
data used. This can be understood from the fact that this 
approach incorporates most physical knowledge. 

In a second set of Monte Carlo simulations, the sen¬ 
sitivity to the initialization of the hyperparameter opti¬ 
mization is analyzed for the three different methods. Only 
one simulated data set is used but the hyperparameter 
optimization is started in 30 randomly selected sets of 
hyperparameters 6q for a varying length of the training 
data. The hyperparameters are assumed to lie around the 
estimates that are obtained using the same optimization 


strategy as above for 8000 data points. Hence, these sets 
of hyperparameters 6^^^^ are known to results in small 
RMSE values as depicted in Figure 4a. The superscript 
‘a’ on is used to explicitly denote that these sets of 
hyperparameters actually differ between the three differ¬ 
ent approaches. For the approach where the magnetic field 
components are modeled using independent CPs, each 
component results in a set of hyperparameters For 
simplicity, in this approach, 6^^^ is chosen to be the mean 
of these three sets of hyperparameters. 

For each of the three approaches, the initial parameters 
Oq are then assumed to deviate from by at most 70% 
as 

^0 = ^tme[l + 0.7U(-l,l)]. (33) 

The Monte Carlo simulation results are depicted in 
Figure 4b. As can be seen, the approach which models 
the three magnetic field components using separate CPs 
suffers most from local minima. Our proposed model 
using a scalar potential still outperforms the other two. 

6.2 Empirical proof-of-concept data 

To illustrate our approach using real data, we have 
performed an experiment where a number of sensors 
have been moved around in a magnetic environment. 
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The sensor board used is shown in Figure 5. We use 
the magnetometer data from an Xsens MTi (Xsens 
Technologies B.V., http://www.xsens.com). Accurate 
position and orientation information is obtained using 
an optical system. These high-accuracy measurements 
were provided through the use of the Vicon real¬ 
time tracking system (Vicon Motion Systems Ltd., UK, 
http://www.vicon.com) courtesy of the UAS Technologies 
Lab, Artificial Intelligence and Integrated Computer 
Systems Division (AIICS) at the Department of Computer 
and Information Science (IDA), Linkoping University, 
Sweden. 

We obtain measurements while sliding the sensor board 
over a configuration of small tables. To ensure sufficient 
excitation, magnets have been placed in an irregular 
pattern underneath these tables as shown in Figure 5. Two 
different data sets have been collected. Both consist of 
approximately three minutes of data sampled at 100 Hz. 
One data set is used for training, while the second is for 
validation. The data of both the training and validation 
data sets are displayed in Figure 6 a. To give an impression 
of the spatial variation of the magnetic field, the magnetic 
field intensity has been visualized through the colors of 
the data. Note that the magnetometer is calibrated such 
that it has a magnitude of one in a local undisturbed 
magnetic field. 

The magnetometer inside the Xsens IMU measures 
the magnetic field at the different locations. The 
optical measurements are used for two purposes. First, 
the positions from the optical system are used as 
known locations in the GP approach. This is a fairly 
reasonable assumption due to the high accuracy of 
the measurements of the optical system. Second, the 
orientations estimated by the optical system are used 
to rotate the magnetometer measurements from the 
magnetometer sensor frame to the lab frame. This rotation 
of the magnetic field measurements is needed for any 
of the GP methods discussed in this work. To use the 
optical and magnetometer data together, they need to be 
time synchronized. This synchronization is done in post¬ 
processing by correlating the angular velocities measured 
by the optical system and by the gyroscope in the Xsens 
IMU. 

We run Algorithm 1 for the training data set. The 
domain U in the GP model is set as Li = I /2 = 0.6 
and I /3 = 0.1. The actual two-dimensional movement 


Optical marker 


|- Xsens IMU Magnets 



Figure 5. Setup of the experiment. Left: Sensor board that 
was used for collecting the data. Right: Magnetic 
environment used in the experiment shown from below 
comprising small magnets to ensure sufficient excitation of 
the magnetic field. 


is performed in a rectangle of 80 cm x 100 cm and is 
hence well within the domain U. The computed scalar 
potential is shown in Figure 6 b. From this the predicted 
magnetic field measurements can be computed which are 
shown in Figure 6 c. For completeness, the intensity of 
these predicted magnetic field measurements are shown in 
Figure 6 d. Although this quantity is only indirectly related 
to the outcome of the GP approach, it is frequently used 
in the remaining sections because of its easy and intuitive 
visualization. 

By predicting the measurements at the locations of 
the validation data set, it is also possible to compute the 
RMSE on the validation data. In Figure 7 we visualize 
the RMSE as a function of the number of basis functions 
used in Algorithm 1. We also compare to the RMSE 
from a full GP approach, that is using the same GP 
prior but without the Hilbert space approximation scheme 
to speed up the inference. To allow for comparison 
with a full GP approach—which suffers from a high 
computational complexity for large data sets—the data 
has been downsampled to 5 Hz. As can be seen, already 
for around m = 1000 basis functions, the quality of the 
estimates from Algorithm 1 approaches that of the full 
GP approach. 

6.3 Mapping the magnetic fieid in a buiiding 

The third experiment was concerned with estimating a 
map of the magnetic environment inside a building by 
only using a smartphone for the data collection. The 
mapped venue is located on the Aalto University campus. 
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(a) Training data (left) and validation data (right) 




(c) The predicted magnetic field. From left to right: and 


(b) Learned scalar potential 


1 

0.5 
0 

(d) Norm of the magnetic field 



Figure 6. Illustration of the magnetic field data and the results from the GP approach from Algorithm 1 for the experiment 
discussed in Section 6.2. In all figures, the y-axis is —0.5, ...,0.5m. The x-axis is —0.4,..., 0.4 m. The units in the field 
surface plots are arbitrary due to normalization. 



Figure 7. RMSE as a function of the number of basis 
functions used in Algorithm 1 as compared to a validation 
data set. The approximative model approaches the full GP 
approach. 

and a floor plan sketch is shown in Figure 8. For practical 
reasons, we limited our interest to the lobby which is 
approximately 600 m^ in size. 


For the measurements, we used an Apple iPhone 4 
and its built-in 9-dof IMU (3-axis AKM AK8975 
magnetometer). All sensors were sampled at 50 Hz, and 
the data was streamed online to a laptop computer for 
processing and storing. The phone was held at waist- 
height and pointed towards the heading direction. 

For reconstructing the walking path and phone 
orientation, we used a pedestrian dead-reckoning (PDR) 
approach developed at Indoor Atlas (Indoor Atlas Ltd., 
Finland, http://www.indooratlas.com), where only 
accelerometer and gyroscope readings were used—the 
path reconstruction thus being fully independent of the 
magnetometer readings. The alignment to the map and 
drift correction were inferred from a set of flxed points 
along the path during acquisition. Two sample paths are 
shown in Figure 8: one of the three training paths with 
similar routes and the validation path. The reconstructed 
paths were visually checked to match the True’ walking 
paths. 

We covered the walkable area in the lobby with three 
walking paths following the same route in each of them 
(see Fig. 8). The reconstructed paths were approximately 
242, 253, and 302 meters long, respectively. The number 
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Figure 8. A training (red) and validation (blue) free-walking path that was used in the experiment. Trajectories were 
collected by a mobile phone, and the magnetometer data was corrected for gravitation direction and heading using the 
inertial sensors in the device. Walking direction markers are shown every 10 meters. The domain boundaries for the 
reduced-rank method are shown by the dashed line. 

of magnetometer data samples acquired along the paths 
were 9868, 10500, and 12335. Prior to each acquisition, 
the phone magnetometer was calibrated by a standard 
spherical calibration approach. The combined size of 
the training data set was n = 32703. For validation, 
we collected a walking path passing through the venue 
(length 54 m, n = 2340). 

We considered a batch interpolation problem of 
creating a magnetic map of the lobby. The map 
was assumed static over time, and we applied Algo¬ 
rithm 1 to the training data with m = 1024 basis 
functions. The optimized hyperparameters were ^ 

575 (/iT)2, ctIe/^Ie 373 (/iT)2, - 1.87m, and 

^noise ~ 5.53 (/iT)^. The coefficients of the inferred lin¬ 
ear (bias) model corresponding to the linear covariance 
function was (—1.095,12.995, —41.119). 

Figure 9 shows the interpolated magnetic field 
magnitude and vector field components. The overall shape 


of the estimate agreed even when the model was trained 
separately with each of the training paths. To most part, 
the strong fluctuations in the magnetic field are located 
near walls or other structures in the building. The strong 
magnetic field in the open area in the lower right part 
of the fioorplan was identified to most likely be due 
to a large supporting structure on the lower fioor-level. 
We also used the model for predicting the measurements 
along the validation part (see Fig. 8). The component¬ 
wise RMSEs were (2.35 /iT, 3.05 /iT, 2.71 /iT) and mean 
absolute errors (1.72 /iT, 2.42 /iT, 2.03 /iT). Considering 
the measurement noise level of the magnetometer (in 
the magnitude of 1 /iT) and the uncertainty in the PDR 
estimate, the solution is surprisingly accurate. 
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(b) Vector field components 


Figure 9. (a) The magnetic field strength (||f ||) interpolated 
by the scalar potential GP model, (b) The separate field 
components of the estimate. 


6.4 Online mapping 

Finally, we demonstrate the power of sequential updating 
and time-dependent magnetic field estimation. The 
mapping was performed by a lightweight and inexpensive 
mobile robot equipped with a magnetometer, and the task 
was to obtain an estimate of the magnetic environment 
of an indoor space by re-calculating the estimate in an 
online fashion. In the second part of the experiment the 
magnetic environment was abruptly changed during the 
experiment, and the aim was to catch this phenomenon by 
spatio-temporal modeling. 

We used a robot for collecting the data. The 
robot was built on a DiddyBorg (PiBorg Inc., UK, 
http://www.piborg.org) robotics board, controlled by a 



Smartphone 


Optical marker 


Invensense IMU 


Trivisio IMU 


DiddyBorg robot board 


Figure 10. The robot was built on a DiddyBorg robotics 
board and controlled by a Raspberry Pi single-board 
computer. The three sensors (Invensense, Trivisio, and 
smartphone) providing magnetometer readings were 
mounted on the top. The reference locations were provided 
by a Vicon optical tracking system. 


Raspberry Pi 2 (model B) single-board computer (Rasp¬ 
berry Pi Foundation, UK, http://www.raspberrypi.org). 
For this example, we controlled the robot over Bluetooth 
with a joystick. 

The robot was equipped with a 9-dof MPU-9150 
Invensense IMU unit that was sampled at 50 Hz. 
The data were collected and stored internally on the 
Raspberry Pi. For additional validation, a Trivisio 
Colibri wireless IMU (TRIVISIO Prototyping GmbH, 
http://www.trivisio.com/), sampled at 100 Hz, and a 
Google Nexus 5 smartphone (AKM AK8963 3-axis 
magnetometer), sampled at 50 Hz, were also mounted 
on the robot for checking the quality of the Invensense 
IMU data. To reduce disturbances caused by the robot, 
the sensors were mounted on an approximately 20 cm 
thick layer of Styrofoam. During post-processing the data, 
the sensor positions and alignments on the robot were 
corrected for. 

High-accuracy location and orientation reference mea¬ 
surements were provided through the use of the Vicon 
(Vicon Motion Systems Ltd., UK, http://www.vicon.com) 






































18 


real-time tracking system courtesy of the UAS Technolo¬ 
gies Lab, Artificial Intelligence and Integrated Computer 
Systems Division (AIICS) at the Department of Computer 
and Information Science (IDA), Linkoping University, 
Sweden. The location measurements could alternatively 
be recovered by odometry and heading information pro¬ 
vided by the robot, but the interest in this experiment was 
rather to focus on the interpolation of the magnetic field, 
not the path estimation. 

The task was to map the magnetic field inside a marked 
region roughly 6 m x 6 m in size. The size of the region 
was limited by the field of view of the Vicon system. The 
magnetometers were calibrated in the beginning of the 
measurement session by rotating the robot around all of its 
axes. A standard spherical calibration approach was used. 
Due to limits in acquisition length of the Vicon system, we 
captured the data in parts, each roughly three minutes in 
length. The magnetic environment remained unchanged 
for the first five data sets (paths shown in Figs, llc-llg), 
and later on changes in the field were initiated by bringing 
in large metallic toolbox shelves. 

In the first part of the experiment, for interpolating 
the magnetic field we used the sequential reduced- 
rank scalar-potential approach presented in Algorithm 2. 
We assumed the magnetic environment stationary, and 
performed sequential updates in an online fashion. For 
practical reasons the calculations were done off-line, but 
the algorithm is fast enough for running in real-time. The 
rank of the approximation was fixed to m = 1024. 

The length-scale, magnitude, and noise variance 
hyperparameters were learned from the first two data sets 
(n = 17980 vector valued observations) by maximizing 
with respect to marginal likelihood (see Alg. 1). 
The obtained values were ^se ~ 0.32 m, ^ 

287 (yuT)^ and ^ 3.27 (/iT)^. The linear model 
magnitude scale parameter was fixed to = 500 (/iT)^. 
The noise model is not only capturing the sensor 
measurement noise, but the entire mismatch between 
the data and the model. This explains the rather large 
noise variance. We also checked, that the hyperparameter 
estimates remained stable when optimized using the rest 
of the data. 

Figure 11 shows the results for the static magnetic 
field experiment. The estimate was updated continuously 
five times a second, and we show five snapshots of the 


evolution of the magnetic field estimate in Figures llc- 
llg (vector field magnitude shown in figures). These 
snapshots also show the path travelled since the previous 
snapshot. The alpha channel acts as a proxy for 
uncertainty; the marginal variance of the estimate is 
giving the degree of transparency. The final magnitude 
estimate—after iterating through all the n = 43029 
observations—is shown in Figure 11a together with the 
vector field components in Figure 11b. 

The frontal part of the mapped region shows strong 
magnetic activity, whereas the parts further back do not 
show as strong fields. Inspection of the venue suggested 
metallic pipelines or structures in the floor to blame (or 
thank) for these features. In this particular case most parts 
of the effect is seen in the x-component. We repeated 
the reconstruction with data collected from the Trivisio 
and smartphone sensors, and the results and conclusions 
remained unchanged. 

The last part of the experiment was dedicated to 
dynamical (time-dependent) modeling of the magnetic 
field. We used all the data from the first part of the 
experiment to train a sequential model and used that as 
the starting point for changing the field (t = 0). During 
acquisition of data while the robot was driving around, 
we brought in two metallic toolbox shelves: first a larger 
toolbox shelf on wheels (Fig. 12c) and then a smaller 
box (Fig. 12d). We acquired altogether some 300 s of data 
(n = 15513) of the changed environment. 

For encoding the assumptions of a changing mag¬ 
netic field, we used Algorithm 2. The additional hyper¬ 
parameter controlling the temporal scale was fixed to 
Aime = 1 hour, thus encoding an assumption of slow local 
changes. This choice is not restrictive, because the data is 
very informative about the abrupt changes. 

Figure 12a shows the evolution of the magnetic field 
components for one fixed location (indicated by a red 
cross in the figures). The two toolboxes induce clear 
changes in the local anomaly field, but the effects are 
restricted to the immediate vicinity of the boxes. Thus 
the spatio-temporal model only gains information about 
the changed field, when the robot passes by the location 
of interest. This effect is clearly visible around t = 20 s, 
and later on around t = 70 s and t = 130 s. After this 
the estimate stabilizes and only drifts around for the 
remaining time. Even though the changes in the field 
components appear clear in Figure 12a, they are only 
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^-component 




(a) Final interpolated magnitude map 


(b) Vector field 



(c) Snapshot 1 (d) Snapshot 2 (e) Snapshot 3 (f) Snapshot 4 (g) Snapshot 5 


Figure 11. The GP interpolation task the robot was faced with. The final interpolation outcome of the magnitude field is 
shown in (a), and the different vector field components are shown in (b). Snapshots along the temporally updating field 
estimate are shown in (c-g) together with the path travelled since the previous update. The marginal variance (uncertainty) 
is visualized by the degree of transparency. 


around 2 /iT and thus only account for a variation of about 
2% in the scale of the entire field visualized in Figure 1 lb. 
Inducing more noticeable changes in the magnetic field 
would require moving around larger structures (say an 
elevator). Yet, even changes this small can be tracked by 
the modeling approach. 

7 Discussion 

Mapping of the magnetic field by Gaussian processes 
appears widely used in the robotics community. This 
paper has aimed both to present a new efficient method 
for mapping, but also provide a study of best practices 
in using GPs in this context. Thus we did go through 
three different approaches for formulating GP priors for 
the magnetic field (independent, shared hyperparameters, 
and curl-free/scalar potential). These three models differ 
in the amount of prior knowledge encoded in the model. 


and the more information available in the prior, the better 
the interpolation and extrapolation capabilities in the 
model—as long as the data agrees with the assumptions. 
This was also demonstrated in Figure 3 and in Section 6.1. 

The methods presented in this paper are related to 
the use of Gaussian random field priors in inverse 
problems (Tarantola 2004; Kaipio and Somersalo 2005). 
This connection has been explored from various points of 
views (cf. Sarkka 2011). However, the machine learning 
(Rasmussen and Williams 2006) way of interpreting 
the GP priors indeed brings something new on the 
table—we are explicitly modeling the uncertainty in the 
field by using a stochastic model, which has interesting 
philosophical implications. The formulation of the prior 
through a GP covariance function provides both an 
intuitive and theoretically justified way of encoding the 
information. 
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Figure 12. (a) Evolution of the magnetic field at one spatial 
location over a time-course of 200 seconds, (b-d) The blue 
metallic toolboxes are brought in at the beginning of the 
experiment. The abrupt changes in the field estimate 
corresponds to time instances when the robot has passed 
the toolboxes and gained information about the changed 
environment. 


Combining models from physics with GPs has also 
been studied under the name Latent force models by 
Alvarez and Lawrence (2009); Alvarez et al. (2013). The 
connection of these models with spatio-temporal Kalman 
filtering was studied, for example, in Hartikainen and 
Sarkka (2011); Hartikainen et al. (2012); Sarkka and 
Hartikainen (2012); Sarkka et al. (2013). However, the 
Kalman filtering approach itself dates back to Curtain and 
Pritchard (1978) in the same way that GPs date back to 
O’Hagan (1978). 

The scalar potential approach is not the only way to 
build the model. It would also be possible to include 
disturbance in the model, which would model the effect 
of free currents in the area or its boundaries. Furthermore, 
more complicated assumptions of the temporal time¬ 
changing behavior of the magnetic field could be included 


in the temporal covariance function. For example, various 
degrees of smoothness or periodicity could be included in 
the framework. 

This paper has been considering the ‘M’ (mapping) 
part in SLAM. The online mapping scheme in this paper 
could be combined by simultaneous localization. As seen 
in the experiments, this appears feasible and provides an 
interesting direction for further research. 

8 Conclusion 

Small variations in the magnetic field can be used as 
inputs in various positioning and tracking applications. 
In this paper, we introduced an effective and practically 
feasible approach for mapping these anomalies. We 
encoded prior knowledge from Maxwell’s equations 
for magnetostatics into a Bayesian non-parametric 
probabilistic model for interpolation and extrapolation of 
the magnetic field. 

The magnetic vector field components were modeled 
jointly by a Gaussian process model, where the prior 
was associated directly with a latent scalar potential 
function. This ensures the field to be curl-free—a justified 
assumption in free spaces. This assumption couples the 
vector field components and additionally encodes the 
assumption of a baseline field with smooth small-scale 
variations. We also presented connections to existing 
formulations for vector-valued Gaussian process models. 

In addition to constructing the model, we also presented 
a novel and computationally efficient inference scheme 
for interpolation and extrapolation using it. We built 
upon a Laplace operator eigenbasis approach, which falls 
natural to the formulation of the model. The inference 
scheme ensures a linear computational complexity with 
respect to the number of observations of the magnetic 
field. We also extended the method to an online approach 
with sequential updating of the estimate and time- 
dependent changes in the magnetic field. 

We presented four experiments demonstrating the 
feasibility and practicality of the methods. A simulated 
experiment showed the benefit of including additional 
knowledge from physics into the model, and a simple 
proof-of-concept example demonstrated the strength of 
the approximation scheme in solving the model. Two real- 
world use cases were also considered: we mapped the 




(c) t = 6 s 


















Solin, Kok, Wahlstrom, Schon, Sarkka 


21 


magnetic field in a building on foot using a smartphone, 
and demonstrated online mapping using a wheeled robot. 
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Appendices 

A Background on the magnetic field 

In this section we go through the magnetostatic 
theory in sufficient detail to transform them into a 
priori assumptions used in the paper. We also present 
the connection between curl-free kernels and scalar 
potentials. 

A. 1 Magnetostatics 

Magnetic fields (see, e.g., Jackson 1999; Vanderlinde 
2004; Griffiths and College 1999) are the effect caused 
by moving electric charges and the magnetic moments 
of elementary particles, and in special relativity they 
are governed by the intercoupling of electromagnetic 
interaction. 

In electrostatics, stationary charges produce an electric 
field E that is constant in time, whereas in magnetostatics 
steady currents produce a constant magnetic field. 
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The Biot-Savart law relates the magnetic field to the 
magnitude, direction, length, and proximity of an electric 
volume current: 


B(x) 


47r 


J(x') X (x-x') 

||X-X'||3 


dx', 


(34) 


V(V • A) - V^A = -V^A. (38) 

Together with Ampere’s law (36) we get one Poisson 
equation for each component in A: 

V^A = -^lQ3. (39) 


where J is the current density and the physical constant fiQ 
is the permeability of free space. The units are such that 
the magnetic field B comes out in newton per ampere¬ 
meter, or tesla (T). The cgs unit gauss is also commonly 
used in physics (1 tesla = 10^ gauss). This formula 
gives the magnetic field at a location x in terms of an 
integral over the current density J(x'). The integration 
is over primed coordinates, but in the following, the 
divergence and curl are taken with respect to the unprimed 
coordinates. 

Taking the divergence of B will now give us Gauss’s 
law for magnetism (see Griffiths and College 1999, for a 
derivation): 

V • B = 0, (35) 

which means that the divergence of a magnetic field is 
zero. 

Applying the curl to Equation (34) gives Ampere’s law: 

V X B = /io J- (36) 

Equations (35) and (36) are Maxwell’s equations for 
magnetostatics. 

A.2 Magnetic vector potential 

Considering the divergence-free magnetic field makes us 
free to consider a vector potential in magnetostatics: 

B = V X A. (37) 

This re-formulation automatically takes care of V • B = 0 
since the divergence of a curl is always zero. With this 
definition, the magnetic vector potential has a built-in 
ambiguity. To the potential A we can add any curl-free 
field without any effect on B. To resolve this ambiguity, 
we also require that V • A = 0 (without loss of generality, 
this is just the so-called Lorenz gauge condition). With 
this assumption we have that 

VxB = Vx(VxA) = 


The vector potential is, however, not as useful as the 
concept of a potential in electrostatics. Even though 
it simplifies the problem formulation in comparison to 
directly working with the Biot-Savart law, it is still 
vector-valued. 

A.3 Magnetic scalar potential 

Instead of the vector potential, we can also consider a 
scalar potential for the magnetic field. To derive such 
a potential we first separate the current density J into 
magnetization current Jm wndfree current Jf as 

J = Jm + Jf • (40) 

The magnetization current is bounded to the material. On 
a microscopic level it can be seen as the current density 
due the orbiting electrons in the material. Since these 
orbits are connected, the magnetization current Jm is 
divergence-free. In a similar manner as in the previous 
section, also the magnetization current Jm has a vector 
potential 

Jm = V X M, (41) 

which is called the magnetization M. By combining (40) 
and (41) with Ampere’s law (36) we have 

Vx(—-M)=Jf, (42) 

VMo J 

=H 

where we denote the difference between the magnetic 
field B and the magnetization M with H. Commonly 
H is called the H-field (some sources refer to H as the 
magnetic field). The assumption of not being inside any 
magnetic material, implies M = 0, and thus B and H are 
proportional. 

In regions where we have no free current Jf = 0, the 
H-field will be curl-free, which is also stated in (1). 
We assume that the ambient magnetic field is induced 
mainly by magnetization in man-built structures in indoor 
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environments. The free-currents are therefore assumed to 
be small or negligible in our domain (Jf ^ 0). The curl- 
free property finally gives that the H-field also has a scalar 
potential 

H = (43) 

which is the key property that is exploited in the model 
developed in this paper. A Gaussian process prior tailored 
for this scalar potential field is presented in Section 4. 
Based on the assumptions discussed above, the scalar 
potential will be a workable description for the ambient 
field that we are interested in. 


which amounts to 
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Thus we get the realation between the scalar potential GP 
magnitude scale and the curl-free kernel magnitude scale 
hyperparameter: 


B Connection between curl-free kernels 
and scalar potentials 

We aim to establish a connection between curl-free 
kernels and the gradient field of the scalar potential 
GP model. Let H G and H' G be two noise-free 
observations of a curl-free magnetic field at x and x', 
respectively. Provided a scalar potential function cp, the 
covariance function of H = — is then determined by 

Cov(H,H')=Keuri(x,x') (44) 

which, provided (p(x) ^ QP(0, /^se(x, x')), is equivalent 
to 
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